Raport

Biblioteki:

library(reshape2)
library(dplyr)
library(tidyr)
library(knitr)
library(ggplot2)
library(plotly)
library(corrplot)
library(shiny)
library(reshape)
library(randomForest)
library(caret)

Wczytanie danych

Zbiór pochodzi z pliku “sledzie.csv”,zwaiera16 atrybutów opisujących 52582 obserwacji. W zbiorze znajdują się 10094 niekompletnych obserwacji. Wartości puste w zbiorze reprezentowane sa za pomocą ‘?’, a część dziesiętna liczb oddzielona jest kropką. Zbiór posiada nagłówek opisujący nazwy atrybutów. By dane zostały wczytane poprawnie, do funcji read.csv przekazano takie argumenty jak: nazwa pliku, forma reprezentacji wartości pustych, że zbiór zawiera nagłówek oraz listę prezentującą klasy danych atrybutów.

Powtarzalność wyników

W celu uzyskania powtarzalności wyników, ustawiono ziarno generatora na:

set.seed(666)

Analiza zbioru:

Na przestrzeni ostatnich lat zauważono stopniowy spadek rozmiaru śledzia oceanicznego wyławianego w Europie. Do analizy zebrano pomiary śledzi i warunków w jakich żyją z ostatnich 60 lat. Dane były pobierane z połowów komercyjnych jednostek. W ramach połowu jednej jednostki losowo wybierano od 50 do 100 sztuk trzyletnich śledzi. Dane dostępne są pod adresem: http://www.cs.put.poznan.pl/alabijak/emd/projekt/sledzie.csv.

Zbiór danych zawiera 16 atrybutów:  

length - długość złowionego śledzia [cm]
cfin1 - dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 1]
cfin2 - dostępność planktonu [zagęszczenie Calanus finmarchicus gat. 2]
chel1 - dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 1]
chel2 - dostępność planktonu [zagęszczenie Calanus helgolandicus gat. 2]
lcop1 - dostępność planktonu [zagęszczenie widłonogów gat. 1]
lcop2 - dostępność planktonu [zagęszczenie widłonogów gat. 2]
fbar - natężenie połowów w regionie [ułamek pozostawionego narybku]
recr - roczny narybek [liczba śledzi]
cumf - łączne roczne natężenie połowów w regionie [ułamek pozostawionego narybku]
totaln - łączna liczba ryb złowionych w ramach połowu [liczba śledzi]
sst - temperatura przy powierzchni wody [°C]
sal - poziom zasolenia wody [Knudsen ppt]
xmonth - miesiąc połowu [numer miesiąca]
nao - oscylacja północnoatlantycka [mb]

Kolejne wiersze reprezentują kolejne obserwacje i są ułożone hronologicznie.

summary(raw_df)

uniq <- raw_df %>% summarise_each(funs(n_distinct(., na.rm = TRUE)))
uniq

Długość śledzia

Długość śledzia jest wyrażona w cm.

plot_length <- ggplot(df,aes(X,length))  + geom_line(alpha=0.3) + stat_smooth(method="auto",size = 2)
plot_length <- plot_length + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) + labs(x="Czas", y="Długość śledzia [cm]") + theme_minimal() 
plot_length
plot_length_hist<- ggplot(df,aes(length)) + geom_histogram(binwidth=0.5, color="white", fill="steelblue") +
  theme_minimal() + labs(x = "Długość", y = "Częstość") + ggtitle("Histogram długości śledzia") 
  
plot_length_hist

Dostępność planktonu - cfin1

plot_plankton1 <- ggplot(df,aes(X,cfin1))  + geom_line(alpha=0.3) + stat_smooth(method="auto",size = 2)
plot_plankton1 <- plot_plankton1 + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) + labs(x="Czas", y="Poziomu dostępności planktonu cfin1") + theme_minimal() 
plot_plankton1
plot_plankton1_hist<- ggplot(df,aes(cfin1)) + geom_histogram(binwidth=0.5, color="white", fill="steelblue") +
  theme_minimal() + labs(x = "Ilość", y = "Częstość") + ggtitle("Histogram dostępności cfin1") 
  
plot_plankton1_hist

Dostępność planktonu - cfin2

plot_plankton2 <- ggplot(df,aes(X,cfin2))  + geom_line(alpha=0.3) + stat_smooth(method="auto",size = 2)
plot_plankton2 <- plot_plankton2 + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) + labs(x="Czas", y="Poziom dostępności planktonu cfin2 od czasu") + theme_minimal() 
plot_plankton2
plot_plankton2_hist<- ggplot(df,aes(cfin2)) + geom_histogram(binwidth=0.5, color="white", fill="steelblue") +
  theme_minimal() + labs(x = "Ilość", y = "Częstość") + ggtitle("Histogram dostępności cfin2") 
  
plot_plankton2_hist

Dostępność planktonu - chel1

plot_plankton3 <- ggplot(df,aes(X,chel1))  + geom_line(alpha=0.3) + stat_smooth(method="auto",size = 2)
plot_plankton3 <- plot_plankton3 + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) + labs(x="Czas", y="Poziom dostępności planktonu chel1 od czasu") + theme_minimal() 
plot_plankton3
plot_plankton3_hist<- ggplot(df,aes(chel1)) + geom_histogram(binwidth=0.5, color="white", fill="steelblue") +
  theme_minimal() + labs(x = "Ilość", y = "Częstość") + ggtitle("Histogram dostępności chel1") 
  
plot_plankton3_hist

Dostępność planktonu - chel2

plot_plankton3 <- ggplot(df,aes(X,chel2))  + geom_line(alpha=0.3) + stat_smooth(method="auto",size = 2)
plot_plankton3 <- plot_plankton3 + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) + labs(x="Czas", y="Poziomu  dostępności planktonu chel2 od czasu") + theme_minimal() 
plot_plankton3
plot_plankton3_hist<- ggplot(df,aes(chel2)) + geom_histogram(binwidth=0.5, color="white", fill="steelblue") +
  theme_minimal() + labs(x = "Ilość", y = "Częstość") + ggtitle("Histogram dostępności chel2") 
  
plot_plankton3_hist

Dostępność planktonu - lcop1

plot_plankton4 <- ggplot(df,aes(X,lcop1))  + geom_line(alpha=0.3) + stat_smooth(method="auto",size = 2)
plot_plankton4 <- plot_plankton4 + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) + labs(x="Czas", y="Poziom dostępności planktonu lcop1 od czasu") + theme_minimal() 
plot_plankton4
plot_plankton4_hist<- ggplot(df,aes(lcop1)) + geom_histogram(binwidth=0.5, color="white", fill="steelblue") +
  theme_minimal() + labs(x = "Ilość", y = "Częstość") + ggtitle("Histogram dostępności lcop1") 
  
plot_plankton4_hist

Dostępność planktonu - lcop2

plot_plankton5 <- ggplot(df,aes(X,lcop2))  + geom_line(alpha=0.3) + stat_smooth(method="auto",size = 2)
plot_plankton5 <- plot_plankton5 + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) + labs(x="Czas", y="Poziom dostępności planktonu lcop2 od czasu") + theme_minimal() 
plot_plankton5
plot_plankton5_hist<- ggplot(df,aes(lcop2)) + geom_histogram(binwidth=0.5, color="white", fill="steelblue") +
  theme_minimal() + labs(x = "Ilość", y = "Częstość") + ggtitle("Histogram dostępności lcop2") 
  
plot_plankton5_hist

Natężenie połowów w regionie - fbar

plot_polow <- ggplot(df,aes(X,fbar))  + geom_line(alpha=0.3) + stat_smooth(method="auto",size = 2)
plot_polow <- plot_polow + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) + labs(x="Czas", y="Natężenie połowów w regionie od czasu - fbar") + theme_minimal() 
plot_polow
plot_polow_hist<- ggplot(df,aes(fbar)) + geom_histogram(binwidth=0.5, color="white", fill="steelblue") +
  theme_minimal() + labs(x = "Natężenie", y = "Częstość") + ggtitle("Histogram natężenie połowu w rejonie -  fbar") 
  
plot_polow_hist

Poziom rocznego narybku - recr

plot_recr <- ggplot(df,aes(X,recr))  + geom_line(alpha=0.3) + stat_smooth(method="auto",size = 2)
plot_recr <- plot_recr + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) + labs(x="Czas", y="Poziom rocznego narybku od czasu - recr") + theme_minimal() 
plot_recr
plot_recr_hist<- ggplot(df,aes(recr)) + geom_histogram(binwidth=160000, color="white", fill="steelblue") +
  theme_minimal() + labs(x = "Roczny narybek", y = "Częstość") + ggtitle("Histogram rocznego narybku - recr") 
  
plot_recr_hist

Roczne natęzenie połowow - cumf

plot_cumf <- ggplot(df,aes(X,cumf))  + geom_line(alpha=0.3) + stat_smooth(method="auto",size = 2)
plot_cumf <- plot_cumf + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) + labs(x="Czas", y="Łączne roczne natężenie połowów w regionie") + theme_minimal() 
plot_cumf
plot_cumf_hist<- ggplot(df,aes(cumf)) + geom_histogram(binwidth=0.5, color="white", fill="steelblue") +
  theme_minimal() + labs(x = "Ilość", y = "Częstość") + ggtitle("Histogram łącznego rocznego natężenia połowów w regionie - cumf") 
  
plot_cumf_hist

Liczba zlowionych ryb - totaln

plot_totaln <- ggplot(df,aes(X,totaln))  + geom_line(alpha=0.3) + stat_smooth(method="auto",size = 2)
plot_totaln <- plot_totaln + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) + labs(x="Czas", y="Liczba złowionych ryb w ramach połowu") + theme_minimal() 
plot_totaln
plot_totaln_hist<- ggplot(df,aes(totaln)) + geom_histogram(binwidth=100000, color="white", fill="steelblue") +
  theme_minimal() + labs(x = "Ilość złowionych ryb", y = "Częstość") + ggtitle("Histogram łącznej liczby złowionych ryb w ramach połowu - totaln") 
  
plot_totaln_hist

Temperatura przy powiechni wody - sst

Atrubut sst zawiera informacje o temperaturze przy powierzchni wody, oraz jest wyrażony w stopniach Celcjusza

plot_sst <- ggplot(df,aes(X,sst))  + geom_line(alpha=0.3) + stat_smooth(method="auto",size = 2)
plot_sst <- plot_sst + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) + labs(x="Czas", y="Temperatura przy powierzchni wody [C]") + theme_minimal() 
plot_sst
plot_sst_hist<- ggplot(df,aes(sst)) + geom_histogram(binwidth=0.5, color="white", fill="steelblue") +
  theme_minimal() + labs(x = "Temperatura", y = "Częstość") + ggtitle("Histogram temperatury przy powierzchni wody - sst") 
  
plot_sst_hist

Poziom zasolenia wody - sal

Atrybut sal opisuje poziom zasolenia wody morskiej. Wartości wyrażone są w Knudsen ppt.

plot_sal <- ggplot(df,aes(X,sal))  + geom_line(alpha=0.3) + stat_smooth(method="auto",size = 2)
plot_sal <- plot_sal + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) + labs(x="Czas", y="Poziom zasolenia wody [ppt]") + theme_minimal() 
plot_sal
plot_sal_hist<- ggplot(df,aes(sal)) + geom_histogram(binwidth=0.05, color="white", fill="steelblue") +
  theme_minimal() + labs(x = "Poziom zasolenia", y = "Częstość") + ggtitle("Histogram poziomu zasolenia wody - sal") 
  
plot_sal_hist

Ostylacja północnoatlantycak - nao

Atrybut nao przedstawia natężenie globalnej cyrkulacji powietrza i wody oceanicznej. Jest zależny od ciśnienia, temperatury, prędkości wiatru i ilości opadów.

plot_nao <- ggplot(df,aes(X,nao))  + geom_line(alpha=0.3) + stat_smooth(method="auto",size = 2)
plot_nao <- plot_nao + theme(axis.text.x=element_blank(),axis.ticks.x=element_blank()) + labs(x="Czas", y="Oscylacja północnoatlantycka") + theme_minimal() 
plot_nao
plot_nao_hist<- ggplot(df,aes(nao)) + geom_histogram(binwidth=0.5, color="white", fill="steelblue") +
  theme_minimal() + labs(x = "Ilość", y = "Częstość") + ggtitle("Histogram oscylacji północnoatlantyckiej - nao") 
  
plot_nao_hist

Analiza brakujących wartości.

Brakujące wartości występują tylko w 7 atrybutach:
- dostępność planktonu -> cfin1, cfin2, chel1, chel2, lcop1, lcop2,
- temperatura przy powierzchni wody -> sst.
Łącznie niepełnych obserwacji mamy: 10094, czyli jest to 19.1966833 % całego zbiioru. Jest to znacząca część zbioru. Jednak zakładamy, że zbiór danych zawiera chronologiczne ułożone dane, dlatego proste operacje na wartościach pustych, takie jak średnia lub mediana z danego atrybutu, mogą całkowicie przekłamać nam rozkład danych. Ze względu na ilość tych danych usunięcie ich może róWnież spowodować różne skutki. Jednak dzięki temu można uznać, że operacje która została dokonana na wartościach pustych, miała jakikolwiek pozytywny lub negatywny na dalsze przetwarzanie.

for (i in 1:ncol(df)){df[is.na(df[,i]),i] <- mean(df[,i],na.rm = T)}

Korelacja między atrybutami

correlation <- cor(df_no_na)
corrplot(correlation,type="upper",tl.col = "black", tl.srt = 45)

Największą korelację można zaobserwować między parami: chel1 - loop1 ,chel2 - loop2, fbar - cumf, cfin2 - lcop2 oraz cumf - totaln. W stosunku do atrybutu ‘length’ najwyższą korelacje można zaobserwować z sst, nao. Nao jest to Oscylacja Północnoatlantycka, która jest zjawiskiem meteorologicznym. Występuje w obszarze północnego atlantyku i ma wpływ na klimat otaczających go kontynentów. Jej działalność związana jest z cyrkulacją powietrza i wody, co by tłumaczyło zależność długości śledzia od tych dwóch atrybutów na raz. length jest także skorelowany z intensywnościa połowów oraz dostępnością planktonu (chel1 - Calanus helgolandicus gat. 1 i lcop1 - widłonogów gat. 1).

Zmiana rozmiaru śledzia w czasie

partition <- createDataPartition(y=df_no_na$length, p=.05, list=FALSE)
dfPartition <- df_no_na[partition, ]
p <- ggplot(dfPartition, aes(x=X, y=length)) + geom_point() + geom_smooth() + theme_bw()
ggplotly(p)

`geom_smooth()` using method = 'gam' and formula 'y ~ s(x, bs = "cs")'

Regresor

Sekcję próbującą stworzyć regresor przewidujący rozmiar śledzia (w tej sekcji należy wykorzystać wiedzę z pozostałych punktów oraz wykonać dodatkowe czynności, które mogą poprawić trafność predykcji); dobór parametrów modelu oraz oszacowanie jego skuteczności powinny zostać wykonane za pomocą techniki podziału zbioru na dane uczące, walidujące i testowe; trafność regresji powinna zostać oszacowana na podstawie miar R2 i RMSE.

Dane zostały podzielone na 2 zbiory: zbiór treningowy - zawierający 80 % danych, oraz zbiór testowy - zawierający 20 %

afterFeatureSelection <- df_no_na %>% select(length, cfin1, cfin2, chel1, chel2, lcop1, lcop2, fbar, recr, cumf, totaln, sst, sal, nao)

inTraining <- createDataPartition(y=afterFeatureSelection$length,p=0.8,list = F)
training <- afterFeatureSelection[inTraining,]
testing <- afterFeatureSelection[-inTraining,]
X_train = select(training, -length)
y_train = training$length
X_test = select(testing, -length)
y_test = testing$length

COS O MODELU

fit <- train(X_train,y_train,
            method = 'rf',
             trControl = ctrl,
             metric = "RMSE",
             tuneGrid=rfGrid,
             importance = TRUE,
             ntree=20)
ggplot(fit) + ggtitle("Optymalizacja parametru mtry na podstawie miary RMSE") + theme_bw()

Sprawdzenie modelu na danych testowych

predictions <- predict(fit, newdata = X_test)
modelValues <- data.frame(obs = y_test, pred = predictions)
fit$results

Model regresyjny został poddany testom na danych testowych. Osiągnięto wynik w postaci miar RMSE i Rsquared:

xs <- 1:length(y_test)
dataf <- data.frame(xs,y_test)
reg_plot_test <- ggplot(dataf,aes(dataf$xs,dataf$y_test)) + 
  geom_point(alpha = 0.5, size = 0.7,color="blue") + 
  geom_smooth(aes(y=dataf$y_test,colour = "trend"), size = 2) +
  coord_cartesian(ylim = c(19, 31)) +
  theme_minimal() + theme(legend.position = "none") +
  labs(x = "Nr przykładu", y = "Długość śledzia ") + 
  ggtitle("Wykres długości śledzia od czasu - prawdziwe dane") 
reg_plot_test

xs <- 1:length(predictions)
dataf <- data.frame(xs,predictions)
reg_plot_test <- ggplot(dataf,aes(dataf$xs,dataf$predictions)) + 
  geom_point(alpha = 0.5, size = 0.7,color="blue") + 
  geom_smooth(aes(y=dataf$predictions,colour = "trend"), size = 2) +
  coord_cartesian(ylim = c(19, 31)) +
  theme_minimal() + theme(legend.position = "none") +
  labs(x = "Nr przykładu", y = "Długość śledzia ") + 
  ggtitle("Wykres długości śledzia od czasu - predykcja") 
reg_plot_test

Analiza ważności atrybutów

importance <- varImp(fit, scale = FALSE)
ggplot(importance)

Z wykresu wynika, że najważniejszym atrybutem jest sst - temperatura przy powierzchni wody. Jednak kolejne 3 najważniejsze atrybuty to poziom 3 różnych rodzajów planktonu - lcop2, cfin2 oraz chel2. Jednak trzeba pamiętać o wysokiej korelacji między cfin2 a loo2, co mogło zadecydować, że oba atrybuty się znalazły jako najważniejsze.

Zależność długości śledzia od temperatury przy powiechni

temp_df <- df_no_na %>% select(X,sst,length)
temp_df2 <- melt(temp_df,id.vars = 'X')
fin_plot_tlen <- ggplot(temp_df2,aes(x = X, y=value)) + 
  geom_point(size = 0.5, alpha = 0.5, color="grey") +
  geom_smooth(size = 2) +
  facet_grid(variable~. ,scales = 'free')  + theme_minimal() + 
  labs(x = "Czas", y = "Temperatura / Długość") + 
  scale_x_continuous() + 
  ggtitle("Wykresy temperatury oraz długości śledzia od czasu") 
fin_plot_tlen

df_planct <- df_no_na %>% select(X,lcop2,cfin2,chel2,length)
df_planct2 <- melt(df_planct,id.vars = 'X')
fin_plot_tlen <- ggplot(df_planct2,aes(x = X, y=value)) + 
  geom_point(size = 0.5, alpha = 0.5, color="grey") +
  geom_smooth(size = 2) +
  facet_grid(variable~. ,scales = 'free')  + theme_minimal() + 
  labs(x = "Czas", y = "Plantkon / Długość") + 
  scale_x_continuous() + 
  ggtitle("Wykresy ilości planktonu oraz długości śledzia od czasu") 
fin_plot_tlen

Wnioski

Interaktywny Wykres

W celu wygenerowania interaktywnego wykresu za pomocą biblioteki shiny, dane z kolumny lenght zostały pogrupowane po 50 kolejnych elementów, a z takich grup policzona średnia. Wykonano to by “suwak” na wykresie nie był za długi. Niestety biblioteka shiny nie daje możliwości na wykorzystanie jej w statycznym pliku raportu, dlatego kod do wykonania takiego wykresu znajduje się w pliku “library.R”.